#include <kuai/misc/Matrix.hpp>
#include <kuai/mol/molalgo.hpp>
#include <kuai/mol/ct.hpp>
#include <kuai/mol/space.hpp>



#ifndef _KUAILIB_SSS_H_
#define _KUAILIB_SSS_H_

namespace kuai {

	template<class MolType1, class MolType2, 
		class PredAtom = std::equal_to<typename MolType1::AtomType>, 
		class PredBond = std::equal_to<typename MolType1::BondType> >
	class Isomophism 
		: public boost::noncopyable
	{
	public:
		typedef typename MolType1::AtomType AtomType;
		typedef typename MolType1::BondType BondType;

	private:
		Isomophism(const MolType1& mol1, const Index* brank1, const MolType2& mol2, const Index* brank2, 
			PredAtom predA = PredAtom(), PredBond predB = PredBond())
			: equalA(predA), equalB(predB)
		{
			assert (mol1.countAtoms() >= mol2.countAtoms());
			assert (mol1.countBonds() >= mol2.countBonds());

			pmolTarget = &mol1;
			pmolPattern = &mol2;

			nAtomT = mol1.countAtoms();
			nAtomP = mol2.countAtoms();
			nBondT = mol1.countBonds();
			nBondP = mol2.countBonds();

			size_t buffersize = nAtomT + nAtomP + nBondT + nBondP;
			if (brank1 == NULL) {
				buffersize += mol1.countBonds();
			}
			if (brank2 == NULL) {
				buffersize += mol2.countBonds();
			}
			buffer.resize(buffersize, 0);

			buffersize = 0;
			if (brank1 == NULL) {
				bond_rankT = &buffer[buffersize];
				RingFinder<MolType1>::getRingRankOfBonds(const_cast<MolType1&>(mol1), const_cast<Index*>(bond_rankT), mol2.countBonds());
				buffersize += mol1.countBonds();
			}
			else {
				bond_rankT = brank1;
			}
			if (brank2 == NULL) {
				bond_rankP = &buffer[buffersize];
				RingFinder<MolType1>::getRingRankOfBonds(const_cast<MolType2&>(mol2), const_cast<Index*>(bond_rankP), mol2.countBonds());
				buffersize += mol2.countBonds();
			}
			else {
				bond_rankP = brank2;
			}

			aamT2P = &buffer[buffersize];
			buffersize += nAtomT;

			aamP2T = &buffer[buffersize];
			buffersize += nAtomP;

			bbmT2P = &buffer[buffersize];
			buffersize += nBondT;

			bbmP2T = &buffer[buffersize];
			buffersize += nBondP;

			nextMapA = nextMapB = 0;

			_aamap.reset(nAtomT, nAtomP, char(-1));
			_bbmap.reset(nBondT, nBondP, char(-1));
		}

		bool search() {
			for (Index i = 0; i < nAtomT; ++i) {
				if (aamap(i, 0)) {
					if (search(i, 0, Index(-1), Index(-1))) {
						return true;
					}
				}
			}
			return false;
		}

		bool search(Index iTarget, Index iPattern, Index iFromT, Index iFromP) {
			Index oldAAIndex = nextMapA, oldBBIndex = nextMapB;

			mapAtom(iTarget, iPattern);

			const AtomType& atomT = pmolTarget->getAtom(iTarget);
			const AtomType& atomP = pmolPattern->getAtom(iPattern);

			Index degreeT = pmolTarget->degree(atomT);
			Index degreeP = pmolPattern->degree(atomP);

			std::vector<Index> indexes(degreeP, 0);
			for (Index i = 0; i < degreeP; ++i) { 
				const AtomType& atomPI = pmolPattern->getNeighborAtom(atomP, i);
				Index iap = pmolPattern->indexOf(atomPI);
				if (iap == iFromP) { 
					continue;
				}
				const BondType& bondPI = pmolPattern->getNeighborBond(atomP, i);
				Index ibp = pmolPattern->indexOf(bondPI);
				if (bbmP2T[ibp] == 0) { 
					Index j = 0;
					for (j = indexes[i]; j < degreeT; ++j) { 
						const AtomType& atomQJ = pmolTarget->getNeighborAtom(atomT, j);
						const BondType& bondQJ = pmolTarget->getNeighborBond(atomT, j);
						Index iaq = pmolTarget->indexOf(atomQJ);
						Index ibq = pmolTarget->indexOf(bondQJ);

						if (aamP2T[iap] == aamT2P[iaq]
							&& bbmP2T[ibp] == bbmT2P[ibq]
							&& aamap(iaq, iap) && bbmap(ibq, ibp))
						{
							mapBond(bondQJ, bondPI);
							if (aamP2T[iap] > 0) 
							{ 
                                break;
							}
							else if (search(iaq, iap, iTarget, iPattern)) {
								break;
							}
							else { 
								restoreBond();
								continue;
							}
						}
					}	// for (j = 0; j < degreeT; ++j)

					if (j == degreeT) { 
						restore(oldAAIndex, oldBBIndex);
						Index k = i-1;
						for (; k < degreeP; --k) { 
							++indexes[k];
							if (indexes[k] < degreeT) { 
								break;
							}
						}
						if (k == Index(-1)) { 
							return false;
						}
						for (k = k+1; 0 < k && k < degreeP; ++k) { 
							indexes[k] = 0;
						}
						i = Index(-1);

						mapAtom(iTarget, iPattern);
					}
				}	// if (bbmP2T[pmolPattern->indexOf(bondPI)]  == 0)
			}	// for (Index i = 0; i < degreeP; ++i)
				
			return true;
		}

		Index mapAtom(Index iQ, Index iP) { 
			assert (aamT2P[iQ] == 0);
			assert (aamP2T[iP] == 0);

            aamT2P[iQ] = aamP2T[iP] = ++nextMapA;
			return nextMapA-1;
		}

		Index mapBond(const BondType& bondQ, const BondType& bondP) { 
			Index iQ = pmolTarget->indexOf(bondQ);
			Index iP = pmolPattern->indexOf(bondP);
			assert (bbmT2P[iQ] == 0);
			assert (bbmP2T[iP] == 0);

			bbmT2P[iQ] = bbmP2T[iP] = ++nextMapB;
			return nextMapB-1;
		}

		void restore(Index oldAAIndex, Index oldBBIndex) { 
			restore(aamT2P, oldAAIndex, nAtomT);
			restore(aamP2T, oldAAIndex, nAtomP);
			restore(bbmT2P, oldBBIndex, nBondT);
			restore(bbmP2T, oldBBIndex, nBondP);
			nextMapA = oldAAIndex;
			nextMapB = oldBBIndex;
		}

		void restore(Index array[], Index oldIndex, Index nLen) { 
			for (Index i = 0; i < nLen; ++i)
			{
				if (array[i] > oldIndex) { 
					array[i] = 0;
				}
			}
		}

		void restoreBond() { 
			nextMapB = nextMapB - 1;
			restore(bbmT2P, nextMapB, nBondT);
			restore(bbmP2T, nextMapB, nBondP);
		}

		bool aamap(size_t i, size_t j) { 
			if (_aamap[i][j] == -1) { 
				const Atom& atomI = pmolTarget->getAtom(i);
				const Atom& atomJ = pmolPattern->getAtom(j);
				_aamap[i][j] = char(0);
				if (equalA(atomI, atomJ)) {
					if (pmolTarget->degree(atomI) >= pmolPattern->degree(atomJ)) {
						_aamap[i][j] = char(1);
					}
				}
			}
			
			return _aamap[i][j] != 0;
		}
		bool bbmap(size_t i, size_t j) { 
			if (_bbmap[i][j] == -1) { 
				const Bond& bondI = pmolTarget->getBond(i);
				const Bond& bondJ = pmolPattern->getBond(j);

				_bbmap[i][j] = char(0);
				if (equalB(bondI, bondJ)) {
					if (bond_rankP[j] == 0 || bond_rankP[j] >= bond_rankT[i]) {
						_bbmap[i][j] = char(1);
					}
				}
			}
			return _bbmap[i][j] != 0;
		}


	private:
		const MolType1* pmolTarget;
		const MolType2* pmolPattern;
		Index nAtomT, nAtomP, nBondT, nBondP;
		const Index* bond_rankT;
		const Index* bond_rankP;
		Index* aamT2P;		// AtomType AtomType map Target to Pattern
		Index* aamP2T;		// AtomType AtomType map Pattern to Target

		Index* bbmT2P;		// BondType BondType map Target to Pattern
		Index* bbmP2T;		// BondType BondType map Pattern to Target

		PredAtom equalA;
		PredBond equalB;

		Index nextMapA, nextMapB;

		std::vector<Index> buffer;

		Matrix<char> _aamap;
		Matrix<char> _bbmap;

	public:

		static bool sss(const MolType1& mol1, const Index* brank1, const MolType2& mol2, const Index* brank2, 
			PredAtom predA = PredAtom(), PredBond predB = PredBond())
		{
			if (mol1.countAtoms() >= mol2.countAtoms() && mol1.countBonds() >= mol2.countBonds()) {
				Isomophism iso(mol1, brank1, mol2, brank2, predA, predB);
				return iso.search();
			}

			return false;
		}


	};

	inline bool sss(const ConnectTable& ctQ, const Index* bondRankQ, const ConnectTable& ctP, const Index* bondRankP) {
		return Isomophism<ConnectTable, ConnectTable>::sss(ctQ, bondRankQ, ctP, bondRankP);
	}
	bool sss(PhaseSpace& ctQ, PhaseSpace& ctP);

}

#endif
